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Abstract 

An analytic leading log formula for the averaged elastic energy loss of a heavy quark is found in the case of a non-static thermalized 
quark gluon plasma. Beyond the mean energy loss the collisional quenching of heavy quarks is considered in an improved pQCD 
transition matrix approach taking into account quantum distribution functions of medium particles, a running coupling and an 
effective screening mass adjusted to HTL calculations. The results for the nuclear modification factor and elliptic flow for RHIC 
t-H energies within this model indicate a large contribution of the collisional energy loss of heavy quarks to their total quenching. 
Thereby hints for a qualitative difference between R AA and were found, such that a linear scaling of the binary cross section 
cannot explain both effects at the same time. Therefore radiative energy loss and more precise initial conditions for the medium 
have to be taken into consideration. 
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1. Introduction 

Jet quenching, as observed in heavy ion collisions at RHIC 
has become an important probe of a deconfined 
state of matter with quasi-free quarks and gluons, the so-called 
OJ quark gluon plasma (QGP). Indeed, this observable allows in- 
| ferences on the properties of the expanding medium. Especially 
jet quenching of heavy flavors charm (c) and bottom (b) could 
be an ideal probe of the QGR Due to the large mas^] of these 
quarks, m » T, they are produced in a very early stage of 
the collision and can later be considered as test particles in the 
(equilibrated) heat bath. 

The quenching of light quarks is well understood as an ef- 
fect of radiative energy loss when propagating through the me- 
dium [5,6]. However, in the case of heavy quarks pure radiative 
£S| energy loss seems not sufficient to explain experimental data 
-H (mainly from electron yields) on heavy flavor quenching [ 7 , 8 , 
[2[T0]|, suggesting surprisingly strong attenuation. There is an 
ongoing debate on the reasons of quenching and energy loss of 
heavy partons in the QGP, which could be possibly clarified by 
^3 analyzing and simulating spectra from heavy flavor decays. An 
additional contribution of suppression, significant in the range 
of small (to medium) momenta, might be the collisional energy 
loss by elastic processes as discussed in (lOl [TTJ [12) and with 
further modifications in (HI [HI US [161 H3. It will turn out 
that in our approach a fairly large collisional energy loss is pre- 
dicted, although some comparisons between collisional and ra- 
diative contributions CUE) indicate another picture where the 
collisional loss should be relatively small (~ 20%) compared 
to the radiative one. These considerations are based on a fixed 
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coupling and use a simple approximation for the screening mass 
in the QGP. In the following work both physical quantities are 
approximated in a more realistic way. Furthermore, references 
|[T3l n~8l point out that some observables might be sensitive to 
the inclusion of binary collisions while the averaged energy loss 
only slightly increases by taking the collisional energy loss of 
parent partons into account, due to different probability distri- 
butions of radiative and collisional processes. 

2. Mean energy loss with QGP flow 

We start with a short overview of the mean energy loss for 
heavy quarks by binary collisions. In a static and thermalized 
plasma case (QED and QCD) a detailed recalculation of the 
energy loss per unit length dE/dx was done in (20l[2l) to revise 
previous calculations by Braaten and Thoma (221I23). The ge- 
neral expression for the mean energy loss of a heavy quark with 
velocity v = p/E and a scattering partner i e {u, d, s, g} from 
the medium is given by 



dEi 
dx 



doj — — oj , 
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as the first moment of the differential transition rate dTi(E)/doj 
with respect to an energy transfer oj = E - E' '. If we are inte- 
rested in computing the constant next to the leading logs of the 
energy loss (see ([4])), we have to decompose the phase space 
of transfered momenta to a soft and a hard part by introducing 
\t*\ as the separation scale, which can be chosen arbitrarily in 
the range of m 2 D <*c \t*\ <*c T 2 , where m 2 D denotes the Debye 
screening mass and T the temperature of the system. There are 
two equivalent ways of computing transition rates either using 
self-energies or matrix elements |22|. That can be seen from 
cutting the one-loop self-energy diagram through the fermion 
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lines. The soft part can be calculated in the so-called hard ther- 
mal loop (HTL) approximation from the self-energies as it was 
done in (20l|2Tl. The calculation of the hard contribution al- 
lows the second approach of computing interaction rates with 
tree-level matrix elements and using a bare propagator instead 
of a dressed one. Here the t-channel long-range gauge inter- 
actions are cut off by imposing a minimum momentum transfer 
\f \ and the expression for the mean energy loss reads: 



dE, 
dx 
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where rii(k) and riffl) in the expressions of the form j£ := 
J d 3 k/(2n) 3 denote the thermal distribution functions. P and 
K stand for the initial and P\ K' for the final four-momenta of 
the collisional partners (P, P' for the heavy quark). The sum 
of the matrix elements \Mi\ 2 runs over all initial and final states 
and is divided by the degeneracy factor d = 6 of the incoming 
heavy quark. In the limit of E » T the equation §2§ can be 
simplified to (2011211 
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here dcri/dt are the Born cross sections l24li for the interaction 
with medium (target) particles of flavor i and degeneracy fac- 
tor d[. A rather long calculation of soft and hard phase space 
contributions leads to the following result for the mean energy 
loss (NLL-formula) of a heavy quark with mass M, which is 
independent of the separation scale \f\ l20ll2Tll : 
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Here c(nf) is a constant beyond the leading logarithms, which 
depends on the number of involved flavors rif = 3 (u,d, s). 

c(n f ) = a-rif + b 0.146 • n f + 0.050 

fl= (f) in(2) 4 + r 

^W(2)-H + H c (5) 
9 / 108 9 

c = — — — y - -1.147 

f(2) ^ 

y ^ 0.577 (Euler's constant) 

This expression of the energy loss is valid for a fixed coupling 
a s . A generalization to the running coupling can be found in 
[21 1. As a remark, the first leading log arises from the t-channel 
interaction and was first anticipated by Bjorken (25). 

In a numerical simulation we use the transition rate (see 
also next chapter) to determine the scattering probability and 
screen the t-channel by a mass term for regulating the infrared 
divergence. This was often done by using the Debye screening 
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Figure 1 : Mean energy loss of a heavy quark as a function of the ingoing mo- 
mentum for two different values of k in comparison to the analytic NLL-formula 





mass, which is not justified as can be seen in fig. [T] where 
the numerical result of the mean energy loss does not converge 
against the analytic formula ([4]). In fact, one has to introduce an 
effective screening mass of the form 



■4n(l + ^W)' 



(6) 



which can contain a fixed or running coupling a(t) and a free 
parameter k. In the case of k = 1 and a fixed coupling ji 2 (t) re- 
duces to the well-known (squared) Debye screening mass m 2 D . 
The procedure of introducing an effective screening mass is le- 
gitimate for the leading-order results because of the previous 
calculations (resumming and renormalizing) on the loop cor- 
rections for gluon exchange processes l26l . Therefore the co- 
efficient k can be calculated by evaluating the interaction rate 
dTi/doj of equation ([T]) with the effective cut-off /j 2 in the t- 
channel and comparing to the result of the NLL-formula ([4}. 
That leads to the following analytic relation for a fixed coupling 
MM 

*=^-0.2. (7) 

2e 

With this value of k we receive a convergence of the numerical 
results (our simulation is described in the next chapter) against 
the analytic formula ([4} as expected, see fig. [T] 

We next derive a relation for the leading log of the mean 
energy loss in the presence of a constant collective QGP flow. 
An interesting case with regard to a heavy ion collision is that 
of a parallel flow in respect to the propagation direction of a 
heavy quark, which interacts with the expanding medium. Here 
we expect a lower energy loss as predicted by the equation ([4]), 
whereas in the antiparallel case the heavy quark should slow 
down faster. As discussed before, the leading logarithmic ac- 
curacy does not require an HTL calculation, which will only 
contribute to a constant next to the leading log. Therefore it is 
sufficient to insert the Debye screening mass squared for the up- 
per bound t max , instead of a decomposition of the phase space. 
The thermal distribution functions rii(k) in the expression of the 
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mean energy loss ([2]) have to be replaced by the Jiittner func- 
tions for bosons and fermions, which depend on the velocity of 
the medium (3 as well as Lorentz-factor y (Lorentz boost of the 
four momentum): 

Hbj = nt\iT ' Hfj = (* A/t ' (8) 

e y{E k -f3-k)/T _ j e y[E k -f3-k)/T + j 

The lower bound on the integration variable t is set by the kine- 
matics and can be approximated for E » M 2 /T in the follow- 
ing way: 
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s = M 2 +2PK+K 2 2Ek 



1 - cos # 



(9(CT)»M 2 , (10) 



where # denotes the angle between p and k. In the same high 
energy limit and \t\ «: s for the leading log region of the t- 
channel, which is dominated by the infrared divergence t — > 0, 
we get 

2M 2 = s + t + u^ s + u, 

( 2\ ( 2 \ (H) 

=> s : = (s - M 2 ) « (M 2 - w) =: -w . 

The evaluation of the t-channel integral with Born cross sec- 
tions l24l and relations ([9]), fTT] ) leads for quarks and gluons at 
leading log accuracy to 
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The u-channel calculation looks similar and requires only a 
change of variables from t to u. Here the dominant contribution 
to the cross section of the Compton scattering is proportional to 
1/u 2 and the relevant kinematic region is given by -s = u min <?c 
u <^ u max = -M 2 s/s as discussed in f20l[2T]|. The remaining 
matrix elements can be neglected, since they contribute only to 
the constant next to the leading log, thus 

— M 2 
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Figure 2: Mean energy loss for a jet of charm quarks pj et = lOGeV as a func- 
tion of the medium velocity J3 according to the NLL formula fTT) for different 
flow directions of the QGP with respect to the momentum vector of a charm 
quark. 



So in the special case of a parallel flow we get 
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here the constant beyond logarithmic accuracy requires the te- 
dious procedure of an HTL calculation but can be approximated 
by the constant c(jif) in the case of a static medium ([5]). 

In the high energy limit E » M, which we use here, the 
term Ey (1 - J3) can be replaced with the Lorentz invariant scalar 
product PiJF. That motivates a simple possibility of calculat- 
ing the general form of dE/dx in the presence of an arbitrary 
but constant flow. Therefore we consider the rest frame (RF) 
of the QGP, where the heavy quark energy is determined by 
the general Lorentz boost of the four momentum vector to this 
frame. Using the equation ([4]), we get immediately the mean 
energy loss in the rest frame of the medium. Now we derive 
a relation for the differential energy loss between RF and a lab 
frame (LF), which is moving with a velocity —p relative to the 
rest frame. The only restriction in our argument is pJF » M at 
the lab frame. We start from the energy loss in RF (primed) and 
express all quantities by corresponding ones in LF (unprimed) 

AE' _ E' 2 -E[ _ j{E 2 -p- p 2 )-y(Ei -P-Pi) 

Af f 2- f i y(h-p-A)-y(h-p-A) 
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_ A£ l~4-^ _ A£ l-y8cos<( J 8,A^|) _ AE 
At 1-/3-^ At l-/3cos< 

(16) 

Here, ± and || denote the transverse as well as parallel contri- 
butions relative to the initial momentum vector p\. The sim- 



pie approximations used here are justified for the mean energy 
loss, since Ap/AE « Ap/Ap and Af/At « Af/Ar are ful- 
filled for /^/^ » M. Furthermore the scattering angle in the 
t- and u-channel for the same energy limit is negligibly small: 
Ap ± /Ap « Ar ± /Ar « 0. Thus the mean energy loss is approxi- 
mately equal in both frames, that yields the leading logarithms 
in the lab frame without any complex calculations. We simply 
replace E of the equation ([4} with the Lorentz invariant quan- 
tity p^ as long as the value of the product is much greater 
as the mass of the heavy quark. The constant next to the lead- 
ing log would require a general evaluation of the integral ^ 
with Jiittner functions and an HTL calculation, which is not be 
investigated here. In general we obtain the following relation 
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Fig. [2] illustrates the mean energy loss according to $FJ\ for 
different flow directions # and a fixed constant c{nf) to repro- 
duce the prediction of the formula ([4} at f3 = 0. We expect from 
the high energy behavior of the mean energy loss (consider ([?]), 
fl6| ) and flT] )) that the constant should be independent of the 
flow velocity p. 

3. MC-simulation with quantum transition rates 

The simulation we use in the following is based on the dou- 
ble differential transition rate T2^p\ p',0) := d 2 T/dp'd6 of bi- 
nary collisions (Jacobian transformation of dTi/daj in ([]])) in a 
thermalized medium of quarks and gluons. The distribution of 
T2 is a function of the momenta p, p' and the scattering angle 
6 (see below). For our purpose T 2 is calculated by a nume- 
rical integration over finite size cells Ap' AO. The divergence 
in the t-channel of the perturbative leading order (Born) cross 
sections is screened with an effective screening mass p 2 (t) as 
described before ([6J. For consistency and accurate predictions 
all a s -dependencies in the screening mass and Born amplitudes 
are calculated with a running coupling as discussed in lfT3lfT4l . 
This procedure generates a two dimensional matrix of transi- 
tion rates. Multiplying the entries of the matrix with a small 
time step At leads to a probability matrix of a Markov chain 
fT3ll . The matrix contains all transition probabilities from an 
ingoing momentum p to another state with momentum p' and 
a scattering angle with respect to the direction of p. A set of 
such probability matrices for different ingoing momenta p form 
the basis of our MC-simulation for the collisional part of heavy 
quark suppression. 

Fig. [3] shows a ratio of numerical to analytical results of 
the mean energy loss (Yf\ with a fixed constant c(rif}. This 
is actually a test of subleading terms, which were neglected in 
the calculation of $T7\ . The figure underlines our prediction 
of the previous chapter that the mean energy loss is dominated 
by the leading log over a long range of j3. Furthermore this 
comparison points out that the beyond leading log formula with 
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Figure 3: Comparison of numerical results for dE/dx to the analytical predic- 
tions for the NLL formula fTT) . The numerical calculations were done within 
the described MC-simulation for a jet of charm quarks with pj et = lOGeV. 



flow produces reasonable results (deviation of - 15%) even at 
rather small energies. 

In our approach the heavy quarks are treated as test particles 
and do not affect the medium. This approximation is justified 
because of the small number of heavy flavor pairs, about 4 com- 
paring to 800 gluons at mid-rapidity 1271 . This allows us to use 
different backgrounds, such as hydro or BAMPS (Boltzmann 
Approach of Multi-Parton Scattering [ 28 ]), for the simulation 
of the local temperature in the medium, which determines the 
collisional transition rate of a thermalized plasma. 

The spacial initialization of heavy quark dijet events is done 
according to the Glauber model (29j [30) and the initial mo- 
mentum distributions are computed with the event generator 
MC@NLO for next-to-leading order calculations [ 3 1 , 32 ]. Here 
the impact of shadowing and Cronin effect is neglected, since 
the modification of heavy quark distributions should be small 
for large momenta pr > 2GeV l33ll . An indirect comparison 
of initial momentum distributions to experimental data can be 
found in ifTTl . 

If at a later stage in our simulation the temperature of the 
medium drops below the critical temperature of T c 200 GeV, 
the heavy quarks pass a fragmentation process to heavy mesons. 
For this we use the Peterson fragmentation | 



N 



(18) 



where N is a normalization constant, z = Pm/Pq the momentum 
fraction (heavy meson to heavy quark) and e c = 0.05, = 
0.005 typical squared mass ratios of involved quarks (light to 
heavy). The following stage of heavy meson decays to electrons 
is calculated with PYTHIA 8.1 l35ll36l . 

One important and free parameter in the described approach 
is the thermalization time r, which has to be fixed by experi- 
mental data or theoretical models. It is not possible to define 
an exact value of the local temperature before equilibrium is 
reached and the calculation of scattering processes of heavy 
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T \T c [GeV] 


To 1 r c [fm/c] 


R[fm] 


dN M /dp} 


set - 1 


0.42 1 0.2 


0.6|5.6 


5.0 


NLO 


set - 2 


0.30|0.2 


1.0|3.4 


6.6 


NLO 



Table 1: Parameter sets for the Bjorken model, R denotes the effective radius of 
the nucleus. 
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Figure 4: Comparison of c + b, D + B (Peterson fragmentation) and electron 
(decay from D + B with PYTHIA) quenching at mid-rapidity in the Bjorken 
model to the observed Raa from heavy flavor decays in the centrality class of 
0% - 10% f7l[38l 



quarks would not be very meaningful. This point will be dis- 
cussed in the next section. 



4. Heavy ion collisions 

For transparency we first consider the particular steps of 
the quenching evolution (for quarks, mesons and electrons) in 
a simple model, namely within the Bjorken model, before we 
investigate more realistic scenarios. Furthermore we general- 
ize the results of fT3lL because these were computed only for 
charm quarks in the transverse plane. Our observable is the 
nuclear modification factor, defined as a ratio of the final trans- 
verse momentum distribution to the distribution in a nucleon- 
nucleon collision, which is scaled with the number of binary 
collision Nb.coii. according to the Glauber assumption 



d 2 N AA /dp t dy 
N bxo ud 2 N NN /dp t dy ' 



(19) 



We replace the power law momentum distributions of lfl"3l by 
the NLO initial momentum distribution and consider the mid- 
rapidity for the same set of parameters (summarized in tab. 
[I]), except the critical freeze-out temperature, which is set to 
T c = 200 MeV. That reduces the life time of the QGP phase 
t c = TqTq/Tc and ensures consistency to later results. Fig. H 
illustrates the combined modification factors of c + b, D + B and 
electron yields for set - 1 and set - 2 (only £-data). From this 
plot it becomes apparent that the c + b spectra are in good agree- 
ment with electron spectra from heavy flavor decays at large p t - 
values (p t >4GeV), but slightly different at intermediate p t and 
underestimate the depletion at smaller p t . Our results for small 



momenta p t < 3 GeV should be considered with some caution, 
as Peterson fragmentation is not suitable here and shadowing 
as well as Cronin effects were neglected. Even in such a sim- 
ple framework we get hints for a significant contribution of the 
collisional energy loss with room for some radiative processes, 
hadronization and further effects. However the radial expansion 
will lead to a shorter life time of the plasma phase, therefore we 
will focus on more realistic dynamical models. 

We can use available partonic or hydro models as input for 
the background medium by extracting the temperature and flow 
data of the discretized space for every time step. We focus 
here on the BAMPS medium, where the effective temperature 
T C eii(x, t) and the flow velocity v ce u(x, t) are extracted locally. 
We then assume in each cell a complete thermodynamical equi- 
librium for the medium. In addition we also use a (3 + 1)- 
dimensional ideal hydrodynamic simulation (SHASTA (37)) as 
a control run to show that the observables are not really sensi- 
tive to differences in these models. Therefore we set the initial 
energy density of the ideal hydro to that of BAMPS by assum- 
ing a thermalization time of r = 0.6 - 0.8 fm/c. These va- 
lues of r are in a good agreement with theoretical predictions 
and experimental data (28| |39) and are similar to the values 
of tab. [T] Fig. [5] confirms that the Raa of charm quarks is 
not sensitive (neglect statistical fluctuations for p t > 7 GeV) 
to a particular starting point of our simulation in the range of 
t = 0.6- 1.0 fm/c. 

Figs. [6] and [7] present our results for the nuclear modifica- 
tion factor in a central and a peripheral collision with an impact 
parameter of b = 8.2 fm, which corresponds to the centrality 
class of 20% - 40%. In both cases we get a quite good agree- 
ment with experimental data, if we keep in mind that the radia- 
tive energy loss is not taken into account. To mimic influence 
of radiative loss an artificial increase of the elastic cross section 
(transition rate) by a factor of K = 3 leads to an overestima- 
tion of the heavy quark quenching. In fact the experimental 
data is enveloped between K = 1 and K = 3 for high transverse 
momenta. Therefore we presume the total cross section of colli- 
sional and radiative processes to vary in this region. That means 
a smaller value of K « 2 in comparison to BAMPS results for 
Raa and v 2 of heavy quarks, which need a K « 4 factor to be 
comparable with experimental data fTTll . We emphasize that 
this observation can be explained from the difference between 
quantum versus Boltzmann distributed particles and from a dif- 
ference of the number of active flavors (BAMPS does not con- 
sider light quarks in reference ifTTIl ). The formula of the mean 
energy loss ^ gets a prefactor of f(2) _1 = 6/n 2 « 0.61 in the 
case of Boltzmann distributed medium particles, which gives a 
rough behavior for high jet energies. 

After the nuclear modification factor R A a we consider now 
the elliptic flow V2, which is given by 



v 2 



P 2 y 



Pt 



(20) 



The elliptic flow of heavy quarks forms at later stages of the 
collision by interactions with the medium. Consequently there 
is only a short and less important time period during the local 
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Figure 5: Sensitivity of Raa in a peripheral nucleus-nucleus collision to diffe- 
rent thermalization times r = 0.6 - 1.0 fm/c, which are used as starting point 
of the simulation. 



thermalization, which cannot be described correctly in our ap- 
proach and is let out (no interactions for heavy quarks). 

Fig. [8] shows V2 in our model in comparison to the expe- 
rimental data. Here we have to emphasize that an artificial in- 
crease of the cross section by a factor of K = 3 can only serve 
as an estimate of the experimental V2 data. Comparing with the 
results of the full BAMPS cascade as published in IT71 . our fac- 
tor for the K- value of 3 is on the same footing to the original 
prediction of BAMPS with K = 4. 

5. Conclusion 

We gave a brief overview of the derivation of the mean en- 
ergy loss for a heavy quark in the QGP and generalized the 
formula to an arbitrary constant medium flow at leading log 
accuracy. These results were numerically tested using our tran- 
sition matrix approach, which includes a running coupling as 
well as a modification of the Debye screening mass inferred by 
HTL considerations and is based on the transition rate for bi- 
nary collisions in a thermalized medium of quarks and gluons. 

In the second part of our work we reconsidered previous re- 
sults on the elliptic flow and nuclear modification factor in the 
same framework of transition matrices. In the case of the nu- 
clear modification factor an increase of the elastic cross section 
by a factor of K < 3 seems to reproduce experimental data of 
electrons from heavy flavor decays, whereas the elliptic flow 
requires K « 3. This discrepancy could be an effect of the ra- 
diative energy loss due to the non-isotropic gluon emission (see 
GUSH) and/or an impact of uncertain initial conditions, which 
could modify the transverse expansion of the medium. 

As a final statement we emphasize that the elastic scattering 
alone cannot explain the sizeable nuclear modification factor 
and elliptic flow of heavy quarks but contributes significantly 
(20% - 50%) to their total quenching. 
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Figure 6: Nuclear modification factor Raa at mid-rapidity \y\ < 0.35 for a cen- 
tral Au + Au collision with a center of mass energy of yfs = 200 GeV (RHIC). 
The curves are calculated within our transition matrix approach for the NLO 
initial momentum distribution of heavy quarks. K denotes an artificial scaling 
factor of the elastic cross section. For comparison, experimental data on elec- 
trons from heavy flavor decays in the centrality class of 0% - 10% are shown 
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Figure 7: Same as fig. [6]but for a peripheral Au + Au collision with an impact 
parameter of b = 8.2 fm. For comparison, experimental data on electrons from 
heavy flavor decays in the centrality class of 20% - 40% are shown [38]. 
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Figure 8: Same as fig.[7]but for elliptic flow. 
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